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Abstract 

An  algorithm  is  presented  which  solves  the  multi- dimensional  ad  vection- diffusion 
equation  on  complex  shapes  to  2^^-order  accuracy  and  is  asymptotically  stable  in  time. 
This  bounded-error  result  is  achieved  by  constructing,  on  a  rectangular  grid,  a  differ¬ 
entiation  matrix  whose  symmetric  part  is  negative  definite.  The  differentiation  matrix 
accounts  for  the  Dirichlet  boundary  condition  by  imposing  penalty  hke  terms. 

Numerical  examples  in  2-D  show  that  the  method  is  effective  even  where  standard 
schemes,  stable  by  traditional  definitions,  fail.  It  gives  accurate,  non  oscillatory  results 
even  when  boundary  layers  are  not  resolved. 
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1  Introduction 


Currently  there  is  a  growing  interest  in  long  time  integration  for  solving  problems  in 
areas  such  as  fluid-mechanics,  aero-acoustics,  electro-magnetics,  material-science,  and  others. 
Clearly,  it  will  be  very  advantageous  if  one  could  formulate  the  spatial  discretization  in  a  way 
which  guarantees  that,  for  the  semi-discrete  formulation,  the  solution-error  norm  is  bounded 
by  the  norm  of  the  truncation  error.  Most,  if  not  all,  existing  algorithms  rely  on  stability 
for  convergence.  However,  even  stable  schemes,  which  at  a  given  time  converge  with  mesh 
refinement  may  have  a  temporally  growing  error,  [1].  This  is  particularly  true  for  hyperbolic 
operators. 

This  paper  considers  2“^-order  accurate  approximations  to  model  linear  advection-diffusion 
equations  in  one  and  more  dimensions,  on  domains  which  may  be  irregular.  By  an  irregular 
domain,  we  mean  a  body  whose  boundary  points  do  not  necessarily  coincide  with  nodes  of 
a  rectangular  mesh. 

In  section  2  we  treat  a  model  “shock-layer”  equation  (linearized  Burger’s  equation), 

tXj  -j-  Q/U^  —  ^  ^  ^5  0  ^  ^  1,  Pi  1. 

K 

We  develop  there  the  theory  for  the  one  dimensional  semi-discrete  system,  resulting  from  the 
spatial  differentiation  used  in  the  finite  difference  algorithm.  Energy  methods  are  used  in 
conjunction  with  “SAT”  type  terms  (see  [1],  [2j),  in  order  to  find  boundary  treatment  and 
“artificial- viscosity-like  terms”,  that  preserve  the  accuracy  of  the  scheme  while  constrain¬ 
ing  an  energy  norm  of  the  error  to  be  temporally  bounded  for  all  t  >  0  by  a  “constant” 
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proportional  to  the  norm  of  the  truncation  error. 

In  section  3  it  is  shown  how  the  methodology  developed  in  section  2  is  used  as  a  building 
block  for  the  multi-dimensional  algorithm,  even  for  irregular  shapes. 

Section  4  presents  numerical  results.  Section  4.1  deals  with  the  steady  state  solution 
to  the  “shock-layer”  equation  for  a  large  range  of  the  “Reynolds  number”,  R.  Oscillations 
that  appear  in  the  numerical  solution  when  using  a  standard  central  finite-differencing,  are 
eliminated  (or  dramatically  reduced)  when  the  bounded-error  algorithm  is  used. 

Section  (4.2)  considers  steady-state  solution  to  a  two  dimensional  scalar  model  to  the 
boundary  layer  equations, 

Ut  +  aUx  +  bUy  =  —Uyy',  R  >>  1 ,  6  <  0, 

IX 

both  for  rectangular  and  trapezoidal  domains.  Again,  the  bounded-error  algorithm  out¬ 
performs  the  standard  scheme  in  ways  described  therein. 

Section  (4.3)  presents  a  time  dependent  example,  modeling  a  boundary-layer  being  ex¬ 
cited  sinosoidially, 

Ut  -1-  aUx  -f  huy  =  ^Uyy  -|-  ufe sin[A;(a:  -  at)]. 

K 

Here,  aside  from  the  usual  performance  criteria,  such  as  error-norms  and  quality  of  the 
velocity  profiles,  we  see  that  the  error-bounded  algorithm  also  has  a  significantly  smaller 
phase  error. 
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2  The  Scalar  One  Dimensional  Case 

Consider  the  scalar  advection-dilFusion  problem 

=  +  +  rL<a;<rR,  t  >  0,  a  >  0,  *  (2.1a) 

u(x,  0)  =  uo(x),  (2.1b) 

=S'L(i),  (2.1c) 

and  f(x,t)  6  C^. 

Let  us  discritize  (2.1)  spatially  on  the  following  uniform  grid: 


Figure  1:  One  dimensional  grid. 

Note  that  the  boundary  points,  x  =  Tz,  and  x  =  Tr,  do  not  necessarily  coincide  with  x^ 
and  XN-  Set  Xj.^i  —  Xj  —  h.,  1  <  j  <  N  —  1]  xx —  Vl  =  0  <  7^  <  1;  Fk  —  xjv  =  7j?^, 

0  <  <  1- 

*The  results  for  the  case  a  <  0  are  found  by  an  analysis  anologus  to  the  one  presented  in  this  section, 
and  are  presented  in  Appendix  I. 
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The  projection  unto  the  above  grid  of  the  exact  solution  u{x^t)  to  (2.1)5  ” 

u{xj^t)  =  u(i).  Let  be  a  matrix  representing  aux  +  ^Uxx:  at  internal  points  without 
specifying  yet  how  it  is  being  constructed.  Then  we  may  write 

4u(t)  =  [i>u(t)  +  B  +  T]  +  m,  (2.2) 

at 

where  T  is  the  truncation  error  due  to  the  numerical  differentiation,  and  fit)  =  f{xj,t), 
1  <  i  <  The  boundary  vector  B  has  entries  whose  values  depend  on  gi,  gR,  'Jl,  Ir 
such  a  way  that  +  B  represents  au^  +  everywhere  to  the  desired  accuracy.  The 

standard  way  of  finding  a  numerical  approximate  solution  to  (2.1)  is  to  omit  T  from  (2.2) 
and  solve 

— v(i)  =  Jbv(t)  +  B  +  f(^)i  (2-3) 

at 

where  v(t)  is  the  numerical  approximation  to  the  projection  u{t).  Subtracting  (2.3)  from 
(2.2)  one  gets  an  equation  for  the  solution  error,  e(t)  =  u{t)  —  v(i), 

4-e  =  De  +  T.  (2.4) 

at 

Our  requirement  for  temporal  stability  is  that  1|  e  ||,  the  L2  norm  of  t.  be  bounded  by 
a  “constant”  proportional  to  /i*”  {rn  being  the  spatial  order  of  accuracy).  Note  that  this 
definition  is  more  severe  than  either  the  G.K.S.  stability  criterion  [3],  or  the  definition  in  [1]. 

It  can  be  shown  that  if  D  is  constructed  in  a  standard  manner,  i.e.,  away  from  the 
boundaries  the  numerical  second  derivative  is  symmetric  and  the  numerical  first  derivative 
is  antisymmetric,  (and  near  the  boundaries  one  uses  “non-symmetric”  differentiation),  then 
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there  are  ranges  of  7^?  and  72,  for  which  D  is  not  negative  definite.  Since  in  the  multi¬ 
dimensional  case  one  may  encounter  all  values  of  0  <  7i,,7jj  <  1,  this  is  unacceptable. 

The  rest  of  this  section  is  devoted  to  the  construction  of  a  scheme  of  2’^'^  order  spatial 
accuracy,  which  is  temporally  stable  for  72,,  72?.  The  basic  idea  is  to  follow  the  procedure  used 
in  [2].  The  present  case  is  more  complicated  due  to  the  difficulty  in  treating  the  advection 
term. 

Note  first  that  the  solution  projection  Uj{t)  satisfies,  besides  (2.2),  the  following  differ¬ 
ential  equation: 

^  =  Du+T,  +  f(t),  (2.5) 

where  now  D  is  indeed  a  differentiation  matrix,  that  does  not  use  the  boundary  values  and 
therefore  Te  T  but  it  too  is  a  truncation  error  due  to  differentiation. 

Next  let  the  semi-discrete  problem  for  v(t)  be,  instead  of  (2.3), 
dv 

—  =  [Dv  -  r2,(ALV  -  gi)  -  Tfi(AfiV  -  g2i)]  -f  f(t),  (2.6) 

where  g2,  =  (1,  •  •  • ,  l)^fl'L(i);  ZR  —  (I5  •  •  •  ?  l)^fl'fi(0?  ^.re  vectors  created  from  the  left  and 
right  boundary  values  as  shown.  The  matrices  A2,  and  Ar  are  defined  by  the  relations: 

A2,u  =  gL  —  T2,,  ArU  =  gR  —  Tr,  (2.7) 

i.e.,  each  row  in  Ai(Ar)  is  composed  of  the  coefficients  extrapolating  u  to  its  boundary  value 
g2,(gH),  at  r2,(r2j)  to  within  the  desired  order  of  accuracy.  (The  error  is  then  T2,(Th)  )• 
The  diagonal  matrices  tl  and  tr  are  given  by 

—  diag  (t2,j , T2,2 5 . . . , T2,j,f ) ,  Tr  diag  {tr^^tr^^  •  •  •  (2.8) 
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Subtracting  (2.6)  from  (2.5)  we  get 


where 


^  =  [De-  tlAlc-  rRARe  +  Tx], 

Cut 


Ti  =  Te  +  TlTl  +  TrTr. 


Taking  the  scalar  product  of  e  with  (2.9)  one  gets: 


(e,  {D  —  tlAl  —  trAr)^  +  (e,  Ti) 
{Z^Me)  +  (c,  Ti). 


(2.10) 


We  notice  that  (e,  Me)  is  (e,  (M  +  M^)e)/2,  where 


M  —  D  —  tjjAi,  —  trAr. 


(2.11) 


If  (M  +  M'^)  can  be  made  negative  definite  then 


(e,(M  +  M^)i)/2< -Co  II  e|p,  (co  >  0). 


(2.12) 


Equation  (2.10)  then  becomes 


€  —Co  II  e  Ip  +(e,Ti) 


and  using  Schwartz’s  inequality  we  get  after  dividing  by  ||  e 


^l|c||<-Co  II  6  II  +  11  Tx 
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and  therefore  (using  the  fact  that  v(0)  =  u(0)) 

II  r||<  II  11^(1  -e-^o*)  (2.13) 

Co 

where  the  “’constanf  ||  Ti  ||m=  max  ||  Ti(r)  ||. 

K  we  indeed  succeed  in  constructing  M  such  that  M  +  is  negative  definite,  with 
Co  >  0  independent  of  the  size  of  the  naatrix  M  as  it  increases,  then  it  follows  from  (2.13) 
that  the  norm  of  the  error  will  be  bounded  for  all  t  by  a  constant  which  is  0{h'^)  where  m 
is  the  spatial  accuracy  of  the  finite  difference  scheme  (2.6).  The  numerical  solution  is  then 
temporally  stable. 

It  can  be  shown  that  as  0,  so  does  co.  When  co  =  0,  the  differential  inequality  is 

I  II  f|l  <  II  T.  II  (2.14) 

leading  to 

||e-i|<||T.  ||„i,  (2.15) 

i.e.,  a  linear  growth  in  time,  a  result  typical  of  hyperbolic  systems.  This  result  can  also  be 
obtained  formally  from  (2.13)  by  letting  co  — 0  for  any  fixed  t. 

The  rest  of  this  section  is  devoted  to  the  task  of  constructing  M  in  the  case  of  m  =  2, 
i.e.,  a  second  order  accurate  finite  difference  algorithm.  We  shall  deal  separately  with  the 
hyperbolic  and  parabolic  parts  of  the  R.H.S.  of  (2.11) 

Let 

M  =  +  o,Mh  =  “  plpAlp  —  trpArp)  +  a{Dii  —  tljjAljj  —  prhArjj).  (2.16) 
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The  parabolic  terms  are  given  by: 


'1-210 
1-210 
0  1-21 
001-21 

=  p  :  (2.17) 

1-2100 
1-2  10 
1  -2  1 
1-21. 

Ti,  =  idiag  0, . . .  O]  =  idiag  o. . . . .  o]  ;  (2.18) 


[0,0, =  ;idiag  [“.O. ■  ■  ■ .  (2  +  1  (^'W) 

^(2  +  7i,)(1 +7l)  -7i(2  +  7n)  2^7l  +  71)  0  ...  0 

;  i  :  i  :  ;  (2.20) 

^(2  +  7l)(1 +71.)  -7i(2  +  7n)  ■^^iL  +  ll)  0  ...  0 

0  ...  0  ^(7i?  +  7n)  -7n(2  +  7n)  +  7n) 

r  .  •  •  • 

0  ...  0  ^(7«  +  7b)  -7fi(2  +  7K)  2^^^  +  1r)0- 1r) 

The  hyperbolic  terms  are  given  by: 


(2.21) 
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Dh  = 


2h 


■  -2  2 

■ 

-1  0 

1 

-1 

0  1 

-1 

0 

1 

0 

-1 

0 

1 

-2 

2  . 

+ 


■  Cl 

■  -1  2  -1 

C2 

0-1  2-1 

C3 

1  -2  0  2  -1 

CjV-2 

1-2  0  2-1 

CjV-1 

1-210 

C/V  . 

1-2  1  . 

where 


and 


”(■  2hc 


■  0 

-1 

1 

> 

1 

-1 

-1  1 

1 

-1 

0  -1 

1 

1  -1 

0 

-1  1 

1 

-1  -1  1 

1  -10. 

> 

(2.22) 


Ck  =  ■  [(cjv  —  ci)k  +  {Nci  —  cjv)], 


(2.23) 


c=-{cx-cn). 


(2.24) 


For  0  >  0  in  (2.1a),  the  left  boundary  is,  for  the  hyperbolic  part,  an  “outflow”  boundary 
on  which  we  do  not  prescribe  a  “hyperbolic  boundary  condition”,  therefore,  in  this  case 
tLh  =  0.  When  a  <  0,  then  trjj  =  0  -  see  the  Appendix  for  details. 
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Here,  with  a  >  0, 


[0,0,...,  ,  4^^] 


(2.25) 


=  0  •  (2-26) 

0  -7r  1  +  7i? 

—'Jr  1  +  7i?  . 

Next  we  shall  show  that  the  parabolic  part  of  M  is  negative  definite.  The  symmetric  part 
of  Mp,  Mp  =  \{Mp  +  Mp),  is  found  using  equations  (2.17)  to  (2.21),  to  be 


37l  —  1 
7L  + 1 

‘2'  — -fL 

2  +  7L 


37z,  —  1  2  —  7l 

7X,  +  1  2  +  7i 


-4  2 


Mp  = 


2  -4  2 


2-4  2 


2  —  7fl 

2  +  7JJ 

37fi  -  1 

7fi  + 1 


2  — 'fR  S'jR  —  1 

2  +  7i?  7i?  +  1 


(2.27) 


We  now  decompose  Mp  as  follows: 


' 

O 

o 

o 

■-4  2 

0  0  0 

2-4  2 

0  0-2  2 

2-4  2 

2-4  2 

Oi 

■*. 

+  (1  -  o:) 

*  .  *  .  *  . 

2-4  2 

2-4  2  0 

2-4  2 

0  2  -2  0  0 

2  -4  . 

0  0  0 

0  0  0. 

-2(1  -  2a) 

®^--)-2a 
IL  +  1 

2-71- 
2  +  71, 

2a 

71-  +  1 

-4(1  -  a) 

2(1  -  a) 

1  —  71- 

2  +  71-' 

2(1  -  a) 

-2(1  -  a) 

0 

0 

0 

0 

-2(1  -  a) 

2(1  -  a) 

2  —  7h 

2  +  7h 

2(1  -  a) 

-4(1  -  a) 

37j?  -  1  ^ 

,  2a: 

7i?  +  1 

2  —  7fl 

2  +  7;? 

2a 

7j?  +  1 

-2(1  -  2a) 

> 

(2.28) 
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We  look  for  1  >  a  >  0  such  that  the  second  and  third  matrices  in  (2.28)  are  non-positive 
definite.  The  first  matrix  in  (2.28)  is  already  negative  definite  by  the  argument  leading  to 
eq.  (2.60),  in  [2].  By  the  same  argument  it  immediately  follows  that  its  largest  eigenvalue  is 
smaller  than  For  0  <  a  <  1,  the  second  matrix  in  (2.28)  is  non-positive  definite,  see 

eq.  (2.63)  k  (2.64)  in  [2].  The  third  matrix  in  (2.28)  has  two  square  3x3  corners  which  are 
negative  for  0  <  a  <  .275.  This  completes  the  proof  that  Mp  is  indeed  negative  definite. 

Next  we  would  like  to  show  that  Mh  =  \{Mh  +  M]j)  is  non-positive  definite.  Using 
equations  (2.22)-(2.26)  we  have 
■  -4  -  2ci  1-1-  2ci  0 

1  -(-  2ci  — 2ci  0 


0  0  0 


0 


0 


2cn  -f  2'^rt^X  -  1  -  2cjv 

-1  -  2cn  -  (1  -1-  7^)'^^^!  +  ^ 


—  (1  + 

i  + 2cm -2(1+ 


We  now  write  Mr  as  the  sum  of  three  “corner-matrices” , 


Mh  =  +  'fR'Hsl, 


where 


(2.29) 


(2.30) 
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— 4  —  2ci  1  “(-  2ci 

0 

1  "I”  2ci  — 2ci 
=  0 


-1  -  (1  +  +  IrTn 


-1  -  (1  +  jr)  4  -  2(1  + 


RRm  =  cn 


2  -2 

-2  2 


(2.31) 


Clearly  is  N.P.D  (non-positive  definite)  for  Vcjv  <  0.  Also,  ttir^  is  N.P.D  for  Ci  >  1/4. 
A  simple  computation  shows  that  is  N.P.D  if  and  tr-  satisfy 


1  +7fi’ 


(<5>0) 


(2.32) 


(/f)  _  1  -  7r{1  —  S) 

“  (1+7KP  ' 


(2.33) 


Thus  we  have  proved  that  Mr  is  indeed  non-positive  definite,  and  therefore  M  =  ^Mp+aMu 
is  negative  definite  for  V4,  a  >  0,  with  its  eigenvalues  bounded  away  below  zero  by  —air^lR, 


0  <  a  <  .275. 
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3  The  Scalar  Two  Dimensional  Case 


We  consider  an  inhomogeneous  advection-diffusion  equation,  with  constant  coefficients, 
in  a  domain  ft.  To  begin  with  we  shall  assume  that  ft  is  convex  and  has  a  boundary  dfl  G  C^. 
The  convexity  restriction  is  for  the  sake  of  simplicity  in  presenting  the  basic  idea;  it  will  be 


removed  later.  The  problem  statement  is: 


du  du  du  d^u  d^u  ,  ,  ^  v,  n 

=  O,-^ - h  b— - 1-  ^1^2  "^^2  0  2  /(^)  y'j  O')  t  >  0,  Ux,  V2  >  0, 

dt  dx  oy  ox^  oy^ 

u(a:,?/,0)  =  uo(a:,y). 


(3.1b) 


u{x,y,'t)\dQ  =  «s(0- 


(3.1c) 


We  shall  refer  to  the  following  grid  representation: 


Figure  2:  Two  dimensional  grid. 


14 


We  have  Mr  rows  and  Mq  columns  inside  Each  row  and  each  column  has  a  discretized 
structure  as  in  the  1-D  case,  see  figure  1.  Let  the  number  of  grid  points  in  the  row  be 
denoted  by  Rk  and  similarly  let  the  number  of  points  in  the  j*'^  column  be  Cj.  Let  the 
solution  projection  be  designated  by  By  U(t)  we  mean,  by  analogy  to  the  1-D  case. 


U(i) 


(^1,1?  ^2,1?  •  •  •  j  '^1,25  ^2,25  •  •  •  7  ^J?2)25  •  •  •  j  '^2,Mjn  •  •  •  i 

(ui,U2,...,umh).  (3.2) 


Thus,  we  have  arranged  the  solution  projection  in  vectors  according  to  rows,  starting  from 
the  bottom  of  fl. 

If  we  arrange  this  array  by  columns  (instead  of  rows)  we  will  have  the  following  structure. 


u(^^(t)  =  (ui,iUi,2,  .  .  .  ,  ;  tt2,l,  «2,2  •  •  .  ,  U2,C2i  UmcA  j  ^<A/c,27  •  •  •  >  ) 


u 


(C) 


>  u 


(C)n 

Me)- 


(3.3) 


Clearly 

U(^)(t)  =  PU,  (3.4) 

where  P  is  an  orthogonal  permutation  matrix,  of  order  i  X  £,  i  being  the  number  of  grid 
points  within  fi. 

The  operator  uid^fdx^  +  a  d/dx  in  (3.1a),  including  the  boundary  terms,  is  represented 
on  the  row  by  whose  structure  is  given  by  (2.16)  and  the  definition  following  it  (see 
(2.17)  through  (2.26)).  Similarly  let  represent  U2d^fdy^  -1-  b  d/dy  on  the  column. 
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With  this  notation,  by  analogy  to  (2.6),  the  two  dimensional  semi-discrete  problem  becomes 


^  +  pTQ^y')  -t-  f(t),  (3.5) 


where  V  is  the  numerical  approximation  of  U; 


■  Mi<** 

1 _ 

xw  = 

;  = 

Mf 

1 

1 _ 

Mil  ■ 

(3.6) 


G(^)  = 


[(^L^gii  +  tr?SRi),  •  •  • ,  +  rjfjgnj, . . . , 

[(4f^gLi + 4f^gj?i)’  •  •  •  ’  ^ 


=  Gf  +  of 


+ 


(3.7) 
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The  subscripts  Bj  (“B”  for  bottom)  play  the  same  role  as  Lk  (“T”  for  left).  The  same 
remark  applied  to  subscripts  Tj  (“T”  for  top)  and  Rk  (“i?”  for  right). 

Note  that  t^(t^)  =  0  when  a  >  0(a  <  0).  Similarly  =  0  when  h  >  0(6  <  0). 

Designating  the  two  dimensional  array  off  errors,  eij,  by  E  =  U  —  V,  the  equation  for  E 
becomes 

^  =  [X(®)  +  P^M^y^P]E  +  T,  (3.8) 

at 

where  T  represents  the  sum  of  the  various  truncation  errors. 

The  time  rate  of  change  of  ||  E  |p  is  given  by 

~  II  E  |p=  (E,  +  P^A^(^)P)E)  +  (E,  T).  (3.9) 

2  at 

By  the  same  argument  that  follow  eq.  (3.15)  in  [2]  it  is  clear  that  the  norm  of  the  error, 
II  E  II,  is  bounded  by  a  constant,  where  the  “constant”  ||  T  j|M=  II  T(r)  ||. 

In  [2],  it  was  shown  that  if  the  domain  Q  is  not  convex  or  simply  connected,  the  above 
results  still  hold.  This  is  also  true  here. 

Note  that  if  ^  ^  =  0  (or  i/i  =  1/2  =  0  in  the  2-D  case)  then  the  differentiation  operator, 

M,  becomes  non-positive  definite.  In  that  case,  it  follows  immediately  from  (3.9)  that  the 
bound  on  the  error-norm  is  not  a  “constant”  but  grows  linearly  in  time. 
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4  Numerical  Examples 

4.1  One  Dimensional  Case 


Here  we  consider  the  problem 


du  1 


t>0,  0  <  a:  <  1, 


(4.1.1) 


u(0,t)  =  1, 
M(l,t)  =  0, 
u(a;,  0)  =  uo(a;). 


The  steady  state  solution  to  (4.1.1)  is: 


,  ,  1  - 


(4.1.2) 


Note  that  R  (=  1  fv)  plays  the  role  of  Reynolds  number  in  this  model  for  a  “linear  shock 
layer” . 

Eq.  (4.1.1)  was  solved  numerically  by  two  methods.  In  one  (referred  to  as  “standard”) 
we  use  central  differencing  for  the  spatial  differentiation,  and  4‘^-order  Runge-Kutta  in  time. 
In  this  “standard”  case,  there  is  no  need  for  special  treatment  at  the  boundaries. 

The  numerical  approximation  v,  in  this  “standard”  case,  satisfies  the  following  finite 
difference  equation: 

=  (0<j<r.)  (4.1.3) 
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with  uo  =  1  and  =  0.  The  solution  to  (4.1.3)  is: 


_  ^  _  2  +  hR 


(4.1.4) 


Notice,  that  if  the  “cell  Reynolds  nunaber,”  Rc  =  hR  >  2,  then  k  <  0  and  the  numerical 
solution,  Vj,  will  be  oscillatory.  If  Rc  <  2  then  we  resolve  the  “shock  layer”  (or  “boundary 
layer”)  and  the  solution  will  be  smooth. 

Numerical  steady-state  solutions  of  (4.1.1)  using  the  “standard  scheme”,  and  using  the 
“bounded-error”  algorithm,  (2,6),  described  above  are  shown  in  figures  3-8  for  Aa:  =  1/100 
and  various  values  of  R.  Both  schemes  were  advanced  to  steady  state  using  4‘^-order  Runge- 
Kutta.  It  is  clear  that  when  Rc  <  2,  both  schemes  give  good  results.  For  Rc  —  10 
{R  =  1000)  both  show  oscillations,  but  the  new  algorithm  approximates  the  exact  solution 
much  better.  When  Rc  =  10^  {R  =  10®),  the  “standard”  numerical  solution  is  useless  while 
the  “bounded-error”  scheme  gives  excellent  results;  in  fact  far  better  than  for  Rc  =  10. 


U 


u 


Figure  3:  Standard  scheme,  Rc  =  2. 


Figure  4:  SAT,  Rc  =  2. 
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Figure  5:  Standard  scheme,  Rc  —  10.  Figure  6:  SAT,  Rc  —  10. 


Figure  7:  Standard  scheme,  Rc  =  Figure  8:  SAT,  Rc  —  1000. 

1000. 

4.2  A  Steady  State  Two  Dimensional  Case 

Here  we  shall  consider  a  steady-state  problem,  which  models,  in  a  way,  the  2-D  boundary 
layer  equations.  The  formulation  is  as  follows:  (The  time  derivative  is  left  in  the  equation, 
since  the  approach  to  steady  state  will  be  via  temporal  advance.) 

Ut  +  au^  +  buy  =  ^Uyy]  t  >  0]  0  <  X  <  V,  0  <  <  1 ,  (4.2.1) 

K 

1  _  pbRy  1  ,  P 

=  I  _  gfeif  +  "  simry,  (4.2.1a) 
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u(a;,0,t)  =  0, 


(4.2.1b) 


u(a:,  l,i)  =  l.  (4.2.1c) 

We  also  take  a  =  1,  and  in  order  to  have  a  growing  “boundary  layer”  on  y  =  0,  we  must  set 
6<0. 

The  analytic  solution  of  this  problem  is: 

.  .  1  —  1  2^  ®  1  •  o 

=  TT^  +  [[ — r 

Figure  9  is  a  3-D  rendition  of  u{x,y)  for  R  =  90,000.  (This  3-D  plot  looks  the  same  to 
the  eye  for  various  —l<b<  — 4/\/R  =  —4/300.)  Figure  10  is  a  plot  of  the  “velocity 
profile”  inside  the  “boundary-layer”  (0  <  y  <  .04)  at  a:  =  .1,.25,.9  and  b  =  —Afy/R.  The 
“bumps”  at  a;  =  .1  and  x  =  .25  may  be  considered  as  “emulating”  results  of  fluid  mechanics 
computation  for  an  incompressible  flow  near  the  entrance  to  a  channel,  see  e.g.  [4]. 

The  numerical  solution  of  (4.2.1)  using  a  standard  central  differencing  scheme  depends 
strongly  on  the  value  of  b  (at  a  given  R).  Figures  11  and  12  show  the  3-D  plot  of  Vj^k  with 
6  =  — 1  and  b  =  —4:fy/R  =  —4/300.  Figs.  13  and  14  show  the  profiles  at  x  =  .1  and  x  =  .9 
for  b  =  —1  and  —355,  respectively.  It  should  be  emphasized  that  the  “peak”  in  figure  11  has 
nothing  to  do  with  the  “bumps”  in  the  exact  solution  (see  figure  10).  The  “peak”  occurs 
way  outside  the  boundary  layer,  and  also  the  amplitude  behavior  with  the  x-coordinate  is 
counter  to  that  of  figure  10.  The  “peak”  is  due  to  a  purely  numerical  oscillation. 

The  same  series  of  plots,  but  as  computed  by  the  new  algorithm,  is  shown  in  figures 
15-18. 
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Figure  17:  SAT,  b=-l.  Figure  18:  SAT,  b  =  -4/300. 

It  should  be  noted  (see  table  1)  that  tbe  “bounded-error”  algorithm  converges  to  steady 
state  (residual  L2  norm  <  10~^^)  an  order  of  magnitude  faster  than  the  standard  scheme 
when  using  the  same  At,  while  cpu-time/iteration  is  about  the  same.  The  standard  scheme 
may  be  run  at  bigger  At  (  by  about  a  factor  of  2)  while  the  SAT  algorithm  was  already  at 
its  maximum  CFL  number.  If  we  let  each  scheme  run  at  its  own  maximum  At  then  the  run 
time  are  about  equal,  but  the  difference  in  errors  remains. 


’time’  to 

L2 

Li  norm 

L2  norm 

Loo  norm 

max  error 

6=  -1 

SAT 

’steady-state’ 

residual 

of  the  error 

of  the  error 

of  the  error 

location 

21.09 

9.911e-14 

8.805e-05 

1.076e-04 

3.108e-04 

45,  46 

Standard 

417 

9.987e-14 

0.485139 

0.674233 

-1.00423 

10,  4 

b  =  -4/300 
SAT 

52.64 

9.943e-14 

1.665e-04 

1.142e-03 

0.01220 

50,  2 

Standard 

416 

9.967e-14 

3.362e-03 

2.447e-02 

-0.2864 

50,  2 

Table  1:  Rectangular  geometry  results. 


We  also  ran  the  same  equations  for  a  non-strictly  rectangular  geometry,  where  the  upper 
boundary  instead  of  being  ?/  =  1  is  t/  =  1  -  (tan^)a:,  where  0  is  the  angle  which  the  upper 
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boundary  makes  with  the  x-axis,  see  figure  19. 


Figure  19:  The  trapezoid  geometry. 

For  many  ^’s  the  results  of  the  performance  of  the  two  schemes  are  unaffected  by  the 
change.  However,  there  are  some  0’s  for  which  the  standard  scheme  converges  to  steady 
state  much  slower  than  before  at  its  own  maximum  allowed  At,  while  the  performance  of 
.the  bounded-error  algorithm  remains  the  same  as  before.  For  example,  see  table  2,  for  the 
case  of  0  =  3.9°.  As  in  [2],  the  point  is  that  for  non-rect angular  geometry  the  distance  that 
a  boundary  is  away  from  a  computational  mode,  'jk,  might  become  extremely  small  and 
this  causes  the  deterioration  in  the  performance  of  the  standard  scheme.  Here  it  is  reflected 
in  the  fact  that  the  standard  scheme  cannot  “support”  the  larger  allowed  At  that  can  be 
achieved  for  the  case  0  =  0.  For  complex  geometries  it  is  very  diflBcult  to  predict  a-priori 
what  range  the  values  of  7  will  take.  The  SAT  methods  (the  bounded  error  algorithm)  is 
insensitive  to  the  variations  in  7  caused  by  the  geometry  of  the  domain. 
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’time’  to  L2  Li  norm  L2  norm  Loo  norm  max  error 

’steady-state’  residual  of  the  error  of  the  error  of  the  error  location 

b  =  -4/300 

SAT  52.56  9.984e-14  1.707e-04  1.156e-03  0.01220  50, 2 

Standard  401.11  9.995e-14  3.448e-03  2.479e-02  -0.2864  50, 2 

Table  2;  Trapezoid  geometry  results. 

4.3  A  2-D  time  dependent  example 


To  check  on  the  temporal  “performance”  of  the  bounded-error  scheme,  we  considered  the 
following  problem: 


ut  +  aU:^  +  buy  =  ^Uyy  +  (Tbsm[k{x-at)];  t>0,  0<x<l,  0<y<l, 

ti 

1  —  bR  bRy  Ib^r?  I  2\  x_  .  .  , 

u(x,  y,  0)  =  j— ^ -t- 2  e  ^  ^  simry +  y(rsmkx, 

,  ,  1  —  bR  bR^  .  •  I  + 

m(0,  y,  t)  =  ^  +  — e  2  sin  xy  -  ya  sm  kat, 


u(a;,0,t)  =  0, 

u{x,l,t)  =  1  4-  cr  siii[A:(x  —  at)]. 


(4.3.1a) 

(4.3.1b) 

(4.3.1c) 
(4.3.1d) 
(4.3. le) 


The  exact  solution  of  (4.3.1)  is: 

.  .  1  —  bR  bRy  (^R?  I  -.2\  X  .  •  XT  f  (A  Q  0^ 

u(x,y,t)  =  - - r^  +  TTT^^  ®  ^iiasmxy-|-yflrsin[A:(a:-at)].  (4.3.2) 

^  1  —  10 

Again  we  take  a  =  1,  i?  =  90,000,  b  —  —l,and  —  The  parameters  <7  and  k  have 
certain  constraints.  If  we  want  u  >  0,  we  must  take  c  <  1.  The  number  of  computational 
nodes,  N,  puts  a  lower  bound  of  2xA  on  the  wave-length,  l/fc,  i.e.,  I  <  k  <  2-kN.  In  the 
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actual  computations  we  used  cr  =  1/2  and  h  =  30.  All  the  plots  for  this  time  dependent  case 
are  shown  for  t  =  10.  Figure  20  shows  a  3-D  plot  of  u{x,  y,  10).  As  in  the  steady-state  case, 
the  plot  looks  the  same  to  the  eye  for  various  — 1  <  6  <  —4fVR  =  —4/300.  Figures  21,  22 
show  the  3-D  plots  of  Vj^k  for  the  standard  and  bounded-error  schemes  respectively.  Figure 

23,  shows  a  x-profile  of  u  at  j/  =  .2,  for  both  schemes  and  the  exact  profile,  for  6  =  1.  Figure 

24,  gives  the  same  profiles  at  j/  =  .8.  These  plots  bring  out  the  differences  in  the  phase  errors 
of  the  numerical  algorithms.  Figures  25-28,  repeat  the  same  information  as  given  in  figure 
21-24,  but  for  b  =  —Ay/R  =  —4/300.  The  efficacy  of  the  bounded-error  algorithm  is  quite 
evident  -  even  when  h  =  —Afy/R,  where  the  norm-errors  away  from  the  boundary  layer  are 
not  dissimilar,  the  phase  error  of  the  right  running  waves  is  quite  a  bit  smaller  in  the  case 
of  the  proposed  present  scheme. 


Figure  20:  Exact  solution.  Figure  21:  Standard  scheme,  b  =  —1. 
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5 


Conclusions 


(i)  A  second  order  method  has  been  developed  which  renders  spatial  second  derivative 
finite  difference  operators  negative  definite.  This  is  not  surprising,  since  negative  defi¬ 
niteness  was  achieved  for  4**^  order  parabolic  operators  in  [2]. 

(ii)  A  second  order  method  has  been  developed  which  renders  spatial  first  derivative  finite 
difference  operators  non-positive  definite.  For  the  case  when  boundary  points  do  not 
coincide  with  grid  nodes  (7  7^  1),  this  is  a  new  result. 

(iii)  The  results  (i)  and  (ii)  allow  us  to  construct  a  solution  operator  for  the  advection 
diffusion  problem  (and,  of  course,  the  diffusion  equation)  which  is  negative  definite, 
thereby  ensuring  asymptotic  temporal  stability. 

(iv)  The  construction  of  these  operators  allows  an  immediate  simple  generalization  to  multi¬ 
dimensional  problems,  on  complex  domains  which  are  covered  by  rectangular  meshes. 
The  proofs  of  the  boundedness  of  the  error-norms  carry  over  rigorously  to  the  (linear) 
multi-dimensional  cases. 

(v)  Numerous  numerical  examples  demonstrate  the  efficacy  of  this  methodology. 
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Appendix  I 

As  in  the  a  >  0  case  the  hyperbolic  terms  are  given  by: 


where 

Ck  =  ^  [(civ  —  ci)k  +  {Nci  —  Cjv)], 


(A.2) 


and 


c=^-{ci-cn). 


(A.3) 
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For  a  <  0  in  (2.1a),  the  right  boundary  is,  for  the  hyperbolic  part,  an  “outflow”  boundary 
on  which  we  do  not  prescribe  a  “hyperbolic  boundary  condition”,  therefore,  in  this  case 

TRjj  =  0, 

and 


TLh  =  7-lf  ^  0, . . . ,  0, 0], 


(A.4) 


1+71/ 

1  +  71- 


-IL 

-7Z-  0 

0 

0 


(A.5) 


Next  we  would  like  to  show  that  Mh  =  \{Mh  +  M]j)  is  non-negative  de&iite, then  aMn 
is  non-positive  definite.  Using  equations  (A.1)-(A.5)  we  have 

-4  -  2ci  -  2(1  +  1  -1-  2ci  -  (H-  'yL)ri^^  +  0 

1  -t-  2ci  -  (1  +  7i-)'r2^^  +  ~2ci  +  0 

0  0  0 


0 


2cn 

— 1  —  Scjv 


(A.6) 


We  now  write  Mr  as  the  sum  of  three  “corner-matrices”, 

Mh  =  +  ^^2  +  ”^^3], 


(A.7) 


0 

-“1  —  2cjv 

4  +  2c;v 
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where 


rriHi  =  Cl 


-2  2 

2-2  0 
0 

0 


mi/2 


-4-2(1+7i,K 


(H) 


1  -  (1  +  0 


1  -  (1  +  +  iLTi 

0 


.(H) 


+2772' 

0 

0 


(H) 


0 

0 


'mHz  — 


0 


2Cjv  — 1  —  2c;v 
—1  —  2cisf  4  +  2civ’ 

Clearly  ttih^  is  N.N.D  (non-negative  definite)  for  Vci  <  0.  Also,  rriHz  is  N.N.D  for 


0 


(A.8) 


cjv  >  — 1/4.  A  simple  computation  shows  that  mjj^  is  N.N.D  if  ti  and  T2  satisfy 


(H) 

r;  ’  = 


2-\-8 
1  +  7l 


(<^  >  0), 


(A.9) 


.(H)  _  1  -  7l(1  -  S) 

(1  +  ihf 


/2  — 


(A.10) 


Thus  we  have  proved  that  Mh  is  indeed  non-negative  definite,  and  therefore  M  =  ^Mp-\-aMif 
is  negative  definite  for  >  0,  with  its  eigenvalues  hounded  away  from  zero  by  — air'll R, 
0  <  a  <  .275,  as  in  the  a  >  0  case  treated  in  the  text. 
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N iim priori  examples  in  2-D  show  that  the  method  is  effective  even  where  standard  scl 
definitions,  fail.  It  gives  accurate,  non  oscillatory  results  even  when  boundary  layers  are 

ion  on  complex  shapes  to 
hiieved  by  constructing,  on 
The  differentiation  matrix 

lemes,  stable  by  tradition 
not  resolved. 

14.  SUBJECT  TERMS 

numerical  solutions;  advection-diffusion  equation; 
a  priori  bounded  errors;  penalty  methods 

15.  NUMBER  OF  PAGES 

36 

16.  PRICE  CODE 

AOS 

17.  SECURITY  CLASSIFICATION 

OF  REPORT 

Unclassified 

18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 

19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

20.  LIMITATION 

OF  ABSTRACT 

NSN  7540-01-280-5500 

Standard  Form  298(Rev.  2-89) 
Prescribed  by  ANSI  Std.  Z39-18 

298-102 


